#include "Noise.h"
#include "wenz_model.h"
#include <iostream>
#include <fstream>
#include <algorithm>
#include <numeric>
#include <cmath>

namespace AlgoEnvNoise {

	std::vector<double> FrequencyRange::generateFrequencies() const {
		std::vector<double> frequencies;
		frequencies.reserve(num_points);

		if (log_scale) {
			double log_min = std::log10(min_freq);
			double log_max = std::log10(max_freq);
			double log_step = (log_max - log_min) / (num_points - 1);

			for (size_t i = 0; i < num_points; ++i) {
				double log_freq = log_min + i * log_step;
				frequencies.push_back(std::pow(10.0, log_freq));
			}
		}
		else {
			double step = (max_freq - min_freq) / (num_points - 1);
			for (size_t i = 0; i < num_points; ++i) {
				frequencies.push_back(min_freq + i * step);
			}
		}

		return frequencies;
	}

	double NoiseSpectrum::getNoiseLevel(double frequency) const {
		if (frequencies.empty() || levels.empty()) {
			return 0.0;
		}

		auto it = std::lower_bound(frequencies.begin(), frequencies.end(), frequency);

		if (it == frequencies.begin()) {
			return levels.front();
		}
		else if (it == frequencies.end()) {
			return levels.back();
		}
		else {
			size_t idx = std::distance(frequencies.begin(), it);
			double f1 = frequencies[idx - 1];
			double f2 = frequencies[idx];
			double l1 = levels[idx - 1];
			double l2 = levels[idx];

			double ratio = (frequency - f1) / (f2 - f1);
			return l1 + ratio * (l2 - l1);
		}
	}

	double NoiseSpectrum::getTotalLevel(const std::pair<double, double>& freq_band) const {
		if (frequencies.empty() || levels.empty()) {
			return 0.0;
		}

		double total_power = 0.0;

		for (size_t i = 0; i < frequencies.size(); ++i) {
			if (frequencies[i] >= freq_band.first && frequencies[i] <= freq_band.second) {
				total_power += std::pow(10.0, levels[i] / 10.0);
			}
		}

		return 10.0 * std::log10(total_power);
	}

	WenzModel::WenzModel(const EnvironmentalParameters& env_params)
		: env_params_(env_params) {
		initializeShippingLevels();
	}

	void WenzModel::setEnvironmentalParameters(const EnvironmentalParameters& env_params) {
		env_params_ = env_params;
	}

	void WenzModel::setShippingLevel(ShippingLevel level, double noise_level) {
		shipping_noise_levels_[level] = noise_level;
	}

	NoiseSpectrum WenzModel::calculateWindNoise(double wind_speed, const FrequencyRange& freq_range) const {
		NoiseSpectrum spectrum;
		spectrum.source_type = NoiseSourceType::WIND;
		spectrum.description = "Wind noise (Wenz model)";
		spectrum.frequencies = freq_range.generateFrequencies();
		spectrum.levels.reserve(spectrum.frequencies.size());

		for (double frequency : spectrum.frequencies) {
			double noise_level = calculateWindNoiseLevel(wind_speed, frequency);
			spectrum.levels.push_back(noise_level);
		}

		return spectrum;
	}

	NoiseSpectrum WenzModel::calculateRainNoise(double rain_rate, const FrequencyRange& freq_range) const {
		NoiseSpectrum spectrum;
		spectrum.source_type = NoiseSourceType::RAIN;
		spectrum.description = "Rain noise (Wenz model)";
		spectrum.frequencies = freq_range.generateFrequencies();
		spectrum.levels.reserve(spectrum.frequencies.size());

		for (double frequency : spectrum.frequencies) {
			double noise_level = calculateRainNoiseLevel(rain_rate, frequency);
			spectrum.levels.push_back(noise_level);
		}

		return spectrum;
	}

	NoiseSpectrum WenzModel::calculateShippingNoise(ShippingLevel shipping_level,
		const FrequencyRange& freq_range) const {
		NoiseSpectrum spectrum;
		spectrum.source_type = NoiseSourceType::SHIPPING;
		spectrum.description = "Shipping noise level " + std::to_string(static_cast<int>(shipping_level));
		spectrum.frequencies = freq_range.generateFrequencies();
		spectrum.levels.reserve(spectrum.frequencies.size());

		for (double frequency : spectrum.frequencies) {
			double noise_level = calculateShippingNoiseLevel(shipping_level, frequency);
			spectrum.levels.push_back(noise_level);
		}

		return spectrum;
	}

	NoiseSpectrum WenzModel::calculateThermalNoise(const FrequencyRange& freq_range) const {
		NoiseSpectrum spectrum;
		spectrum.source_type = NoiseSourceType::THERMAL;
		spectrum.description = "Thermal noise";
		spectrum.frequencies = freq_range.generateFrequencies();
		spectrum.levels.reserve(spectrum.frequencies.size());

		for (double frequency : spectrum.frequencies) {
			double noise_level = calculateThermalNoiseLevel(frequency);
			spectrum.levels.push_back(noise_level);
		}

		return spectrum;
	}

	NoiseSpectrum WenzModel::calculateTotalNoise(double wind_speed, double rain_rate,
		ShippingLevel shipping_level,
		const FrequencyRange& freq_range) const {
		auto wind_noise = calculateWindNoise(wind_speed, freq_range);
		auto rain_noise = calculateRainNoise(rain_rate, freq_range);
		auto shipping_noise = calculateShippingNoise(shipping_level, freq_range);
		auto thermal_noise = calculateThermalNoise(freq_range);

		std::vector<NoiseSpectrum> all_spectra = { wind_noise, rain_noise, shipping_noise, thermal_noise };
		NoiseSpectrum total_spectrum = combineSpectra(all_spectra);
		total_spectrum.description = "Total environmental noise";

		return total_spectrum;
	}

	NoiseSpectrum WenzModel::calculateEnvironmentalNoise(const EnvData::SpatioTemporalCoord & coord,
		ShippingLevel shipping_level,
		const FrequencyRange & freq_range) const {
		auto& data_manager = EnvData::EnvironmentalDataManager::getInstance();

		double wind_speed = 0.0;
		double rain_rate = 0.0;

		try {
			wind_speed = data_manager.getWindData().getValue(EnvData::DataType::WIND_SPEED, coord);
		}
		catch (const std::exception& e) {
			std::cerr << "get wind data fail" << e.what() << "default 5.0 m/s" << std::endl;
			wind_speed = 5.0;
		}

		try {
			rain_rate = data_manager.getRainData().getValue(EnvData::DataType::RAIN_RATE, coord);
		}
		catch (const std::exception& e) {
			std::cerr << "get rain data fail" << e.what() << "default 0.0 mm/h" << std::endl;
			rain_rate = 0.0;
		}

		return calculateTotalNoise(wind_speed, rain_rate, shipping_level, freq_range);
	}

	NoiseSpectrum WenzModel::combineSpectra(const std::vector<NoiseSpectrum>& spectra) {
		if (spectra.empty()) {
			return NoiseSpectrum();
		}

		NoiseSpectrum combined;
		combined.frequencies = spectra[0].frequencies;
		combined.levels.resize(combined.frequencies.size(), 0.0);
		combined.description = "Combined noise spectrum";

		for (size_t freq_idx = 0; freq_idx < combined.frequencies.size(); ++freq_idx) {
			double total_power = 0.0;

			for (const auto& spectrum : spectra) {
				if (freq_idx < spectrum.levels.size()) {
					total_power += std::pow(10.0, spectrum.levels[freq_idx] / 10.0);
				}
			}

			combined.levels[freq_idx] = 10.0 * std::log10(total_power);
		}

		return combined;
	}

	NoiseSpectrum WenzModel::interpolateSpectrum(const NoiseSpectrum& spectrum,
		const std::vector<double>& target_frequencies) {
		NoiseSpectrum interpolated;
		interpolated.source_type = spectrum.source_type;
		interpolated.description = spectrum.description + " (interpolated)";
		interpolated.frequencies = target_frequencies;
		interpolated.levels.reserve(target_frequencies.size());

		for (double freq : target_frequencies) {
			interpolated.levels.push_back(spectrum.getNoiseLevel(freq));
		}

		return interpolated;
	}

	double WenzModel::calculateWindNoiseLevel(double wind_speed, double frequency) const {
		// N_w(f) = 50 + 7.5 * W^0.5 - 20 * log10(f) + correction_factors

		if (wind_speed <= 0.0 || frequency <= 0.0) {
			return -200.0;
		}

		double base_level = WIND_COEFFICIENT + 7.5 * std::sqrt(wind_speed);
		double freq_term = -20.0 * std::log10(frequency);
		double freq_correction = getFrequencyCorrection(frequency);
		double depth_correction = getDepthCorrection(env_params_.depth, frequency);

		return base_level + freq_term + freq_correction + depth_correction;
	}

	double WenzModel::calculateRainNoiseLevel(double rain_rate, double frequency) const {
		// N_r(f) = 35 + 10 * log10(R) + frequency_dependent_term

		if (rain_rate <= 0.0) {
			return -200.0;
		}

		if (frequency <= 0.0) {
			return -200.0;
		}

		double base_level = RAIN_COEFFICIENT + 10.0 * std::log10(rain_rate);

		double freq_factor = 0.0;
		if (frequency > 1000.0) {
			freq_factor = 10.0 * std::log10(frequency / 1000.0);
		}
		else {
			freq_factor = -30.0 * std::log10(frequency / 1000.0);
		}

		return base_level + freq_factor;
	}

	double WenzModel::calculateShippingNoiseLevel(ShippingLevel shipping_level, double frequency) const {
		auto it = shipping_noise_levels_.find(shipping_level);
		if (it == shipping_noise_levels_.end()) {
			return -200.0;
		}

		double base_level = it->second;

		if (frequency <= 0.0) {
			return base_level;
		}

		double freq_correction = 0.0;
		if (frequency < 100.0) {
			freq_correction = 5.0 * std::log10(100.0 / frequency);
		}
		else if (frequency > 1000.0) {
			freq_correction = -10.0 * std::log10(frequency / 1000.0);
		}

		return base_level + freq_correction;
	}

	double WenzModel::calculateThermalNoiseLevel(double frequency) const {
		// N_t(f) = -15 + 20 * log10(f)

		if (frequency <= 0.0) {
			return -200.0;
		}

		return THERMAL_COEFFICIENT + 20.0 * std::log10(frequency);
	}

	double WenzModel::getFrequencyCorrection(double frequency) const {
		if (frequency < 10.0) {
			return -10.0 * std::log10(frequency / 10.0);
		}
		else if (frequency > 10000.0) {
			return -5.0 * std::log10(frequency / 10000.0);
		}
		return 0.0;
	}

	double WenzModel::getDepthCorrection(double depth, double frequency) const {
		if (depth < 10.0) {
			return 3.0; // 
		}
		else if (depth > 1000.0) {
			return -2.0; // 
		}

		// 
		double depth_factor = 1.0 - frequency / 10000.0;
		depth_factor = std::max(0.1, depth_factor);

		return depth_factor * std::log10(depth / 100.0);
	}

	void WenzModel::initializeShippingLevels() {
		shipping_noise_levels_[ShippingLevel::LEVEL_0] = 0.0;
		shipping_noise_levels_[ShippingLevel::LEVEL_1] = 60.0;
		shipping_noise_levels_[ShippingLevel::LEVEL_2] = 70.0;  // 
		shipping_noise_levels_[ShippingLevel::LEVEL_3] = 80.0;  //
		shipping_noise_levels_[ShippingLevel::LEVEL_4] = 90.0;  // 
	}

	// EnvironmentalNoiseCalculator
	EnvironmentalNoiseCalculator::EnvironmentalNoiseCalculator(const EnvironmentalParameters& env_params)
		: wenz_model_(env_params) {
		data_manager_ = std::shared_ptr<EnvData::EnvironmentalDataManager>(
			&EnvData::EnvironmentalDataManager::getInstance(),
			[](EnvData::EnvironmentalDataManager*) {} //
		);
	}

	std::vector<NoiseSpectrum> EnvironmentalNoiseCalculator::calculateRegionalNoise(
		const std::pair<double, double>& time_range,
		const std::pair<double, double>& lat_range,
		const std::pair<double, double>& lon_range,
		ShippingLevel shipping_level,
		const FrequencyRange& freq_range) const {

		std::vector<NoiseSpectrum> regional_spectra;

		// 
		try {
			auto wind_data = data_manager_->getWindData().getRegionData(
				EnvData::DataType::WIND_SPEED, time_range, lat_range, lon_range);
			auto rain_data = data_manager_->getRainData().getRegionData(
				EnvData::DataType::RAIN_RATE, time_range, lat_range, lon_range);

			// 
			size_t data_size = std::min(wind_data.size(), rain_data.size());
			regional_spectra.reserve(data_size);

			for (size_t i = 0; i < data_size; ++i) {
				double wind_speed = wind_data[i];
				double rain_rate = rain_data[i];

				NoiseSpectrum spectrum = wenz_model_.calculateTotalNoise(
					wind_speed, rain_rate, shipping_level, freq_range);

				regional_spectra.push_back(std::move(spectrum));
			}

		}
		catch (const std::exception& e) {
			std::cerr << "Regional Noise Calculation Failed: " << e.what() << std::endl;
		}

		return regional_spectra;
	}

	std::vector<NoiseSpectrum> EnvironmentalNoiseCalculator::calculateTimeSeriesNoise(
		double latitude, double longitude,
		const std::vector<double>& time_points,
		ShippingLevel shipping_level,
		const FrequencyRange& freq_range) const {

		std::vector<NoiseSpectrum> time_series_spectra;
		time_series_spectra.reserve(time_points.size());

		for (double time : time_points) {
			EnvData::SpatioTemporalCoord coord(latitude, longitude, time);

			NoiseSpectrum spectrum = wenz_model_.calculateEnvironmentalNoise(
				coord, shipping_level, freq_range);

			time_series_spectra.push_back(std::move(spectrum));
		}

		return time_series_spectra;
	}

	bool EnvironmentalNoiseCalculator::exportResults(const std::string& filename,
		const std::vector<NoiseSpectrum>& spectra) const {
		std::ofstream file(filename);
		if (!file.is_open()) {
			std::cerr << "Unable to create output file: " << filename << std::endl;
			return false;
		}

		if (spectra.empty()) {
			std::cerr << "No data can be exported" << std::endl;
			return false;
		}

		// 
		file << "spectrum_id,frequency,noise_level,description" << std::endl;

		// 
		for (size_t spec_idx = 0; spec_idx < spectra.size(); ++spec_idx) {
			const auto& spectrum = spectra[spec_idx];

			for (size_t freq_idx = 0; freq_idx < spectrum.frequencies.size(); ++freq_idx) {
				file << spec_idx << ","
					<< spectrum.frequencies[freq_idx] << ","
					<< spectrum.levels[freq_idx] << ","
					<< spectrum.description << std::endl;
			}
		}

		std::cout << "The result has been successfully exported: " << filename << std::endl;
		return true;
	}



	std::vector<std::vector<int>> getShipLevels(const ShipData& shipdata) {
		std::vector<std::vector<int>> ship_levels(shipdata.lat.size(), std::vector<int>(shipdata.lon.size()));

		for (size_t i = 0; i < ship_levels.size(); i++) {
			for (size_t j = 0; j < ship_levels[0].size(); j++) {
				int shipdensity = shipdata.ship_densitis[i][j];

				if (shipdensity > 50000000) {
					ship_levels[i][j] = 4;
				}
				else if (shipdensity > 35000000) {
					ship_levels[i][j] = 3;
				}
				else if (shipdensity > 20000000) {
					ship_levels[i][j] = 2;
				}
				else if (shipdensity > 10000000) {
					ship_levels[i][j] = 1;
				}
				else
				{
					ship_levels[i][j] = 0;
				}
			}
		}

		return ship_levels;
	}

	std::vector<std::vector<double>> getRainFall(const int& rainlevel, const size_t& latSize, const size_t& lonSize) {
		std::vector<std::vector<double>> rain_rates(latSize, std::vector<double>(lonSize));  //(lat, lon)
		double rate;

		if (rainlevel == 0) {
			rate = 0.0;
		}
		else if (rainlevel == 1) {
			rate = 0.8;  // [0.1, 1.5]
		}
		else if (rainlevel == 2) {
			rate = 3.0;  // [1.6, 6.0]
		}
		else if (rainlevel == 3) {
			rate = 13.0;  // [6.1, 20.0]
		}
		else
		{
			rate = 25.0;
		}

		for (auto& lonData : rain_rates) {
			for (auto& rain_rate : lonData) {
				rain_rate = rate;
			}
		}

		return rain_rates;
	}

	TotalNoiseData Noise::AreaTotalNoiseCalculation(const WindData& winddata, const int& rainlevel, const ShipData& shipdata, const NoiseConfig& noiseCfg) {
		std::cout << "\n=== Area Total Noise Calculation ===" << std::endl;

		size_t latSize = winddata.lat.size();
		size_t lonSize = winddata.lon.size();
		auto rain_rates = getRainFall(rainlevel, latSize, lonSize);

		auto ship_levels = getShipLevels(shipdata);

		FrequencyRange freq_range(noiseCfg.min_f, noiseCfg.max_f, noiseCfg.num, noiseCfg.log);

		TotalNoiseData data;
		data.lat = winddata.lat;
		data.lon = winddata.lon;
		data.time = winddata.time;
		data.total_noise = std::vector<std::vector<NoiseSpectrum>>(latSize,
			std::vector<NoiseSpectrum>(lonSize));

		WenzModel wenz;

		for (size_t i = 0; i < latSize; i++) {
			for (size_t j = 0; j < lonSize; j++) {
				double wind_speed = winddata.wind_speeds[i][j];  // m/s
				double rain_rate = rain_rates[i][j];    // mm/h
				int ship_level = ship_levels[i][j];
				ShippingLevel shipping_level = static_cast<ShippingLevel>(ship_level);
				data.total_noise[i][j] = wenz.calculateTotalNoise(wind_speed, rain_rate, shipping_level, freq_range);
			}
		}

		return data;
	}

}
